rm(list = ls())
options(width=110, device = 'quartz')

###################
## Assemble the Data 
####################
library(foreign)
library(plyr); library(dplyr)
library(coin)
library(rdrobust)
library(rms)
library(rdd)
library(ggplot2)
library(reshape2)
library(xtable)
library(tidyr)

myRecode <- function(var, recode.ls, as.factor.result=TRUE, ...) {
  options(useFancyQuotes = FALSE)
  recode.ls <- sapply(recode.ls, function(x) paste(sQuote(x), collapse='='))
  recodes <- paste(recode.ls, collapse=';')
  print(recodes)
  revar <- Recode(var, recodes, as.factor.result=TRUE, ...)
  options(useFancyQuotes = TRUE)
  revar
}
rols <- function (..., data, subset=TRUE, cluster=1:nrow(data),
                  method="huber") {
  mod <- ols(..., data=data[subset, ], x=TRUE)
  robcov(mod, cluster=cluster[subset], method=method)
}
PasteEval <- function(..., print = FALSE) {
  txt <- paste(..., sep = "")
  if (print) { cat("\n", txt, "\n") }
  eval(parse(text = txt), envir = parent.frame())
}
FileName <- function (path=NULL, name=NULL, ext=NULL, replace=FALSE) {
  p <- paste(path, collapse = "")
  n <- paste(name, collapse = "")
  e <- paste(".", ext, sep = "")
  d <- format(Sys.Date(), "%y%m%d") 
  for (l in seq_along(letters)) {
    if (l > 1) old_file_name <- file_name
    file_name <- paste0(p, d, paste(n, letters[l], sep="-"), e)
    if (!any(grepl(file_name, list.files()))) {
      if (replace) file_name <- ifelse(l > 1, old_file_name, file_name)
      break
    }
  }
  cat("\nFile name:", file_name, "\n\n")
  return(file_name)
}

mysw <- function (...) {
  setwd(paste(..., sep = "/"))
}
mypdf <- function(name, ...) {
  date <- format(Sys.Date(), "%y%m%d")
  pdf(paste(paste(name, collapse=""), date, ".pdf", sep=""), ...)
}
tri <- function (x, h, c=0) pmax(0, 1 - abs((x - c) / h))




### Multidimensional RDD Data Prep

### Run data prep code first

load(file= "party_control_data_161121.RData")

round2 = function(x, n) {
  posneg = sign(x)
  z = abs(x)*10^n
  z = z + 0.5
  z = trunc(z)
  z = z/10^n
  z*posneg
}

data$hs_elections_this_year[data$abb=="NJ" & data$year==2014]<-1
data$hs_elections_this_year[data$abb=="VA" & data$year==2014]<-1
stateleg_elections <- data[data$hs_elections_this_year==1,]

stateleg_elections$hs_tot_in_sess[is.na(stateleg_elections$hs_tot_in_sess)]<-stateleg_elections$hs_dem_in_sess[is.na(stateleg_elections$hs_tot_in_sess)] + stateleg_elections$hs_rep_in_sess[is.na(stateleg_elections$hs_tot_in_sess)]

stateleg_elections$margin_abs <- abs((stateleg_elections$hs_dem_in_sess - stateleg_elections$hs_rep_in_sess) /2)
stateleg_elections$margin_abs2 <- round(abs(stateleg_elections$hs_dem_in_sess-(stateleg_elections$hs_tot_in_sess/2)))
stateleg_elections$margin <- ((stateleg_elections$hs_dem_in_sess - stateleg_elections$hs_rep_in_sess) /2)
stateleg_elections$margin2 <- stateleg_elections$hs_dem_in_sess-(stateleg_elections$hs_tot_in_sess/2)
stateleg_elections$margin<-stateleg_elections$margin2
stateleg_elections <- subset(stateleg_elections, election_year>1967)
stateleg_elections$even.sh<-NA
#stateleg_elections$even.sh<- (stateleg_elections$hs_dem_in_sess - stateleg_elections$hs_rep_in_sess) %% 2 == 0 
stateleg_elections$even.sh<- (stateleg_elections$hs_tot_in_sess/2) %% 2 == 0 
stateleg_elections$multimember[stateleg_elections$abb=="VA" & stateleg_elections$year==2014]<-0
stateleg_elections$multimember[stateleg_elections$abb=="NJ" & stateleg_elections$year==2014]<-1
stateleg_elections$multimember[stateleg_elections$abb=="VA" & stateleg_elections$year==2013]<-0
stateleg_elections$multimember[stateleg_elections$abb=="NJ" & stateleg_elections$year==2013]<-1

stateleg_elections$margin_abs2[stateleg_elections$even.sh==T]<-stateleg_elections$margin_abs2[stateleg_elections$even.sh==T]+1

data2<-data
load(file="SLERs1967to2010_2012_05_26.RData")
elections1<-data
data<-data2
elections2<-read.dta(file="SLERs2011to2012_only_2013_05_14.dta")
elections3<-read.dta(file="SLERs2013.dta")
## get the demoratic and republican votes in each election
elections2<-spread(elections2, v21, v23)
elections2$v30_dems<-elections2$'100'
elections2$v31_reps<-elections2$'200'
elections3<-spread(elections3, v21, v23)
elections3$v30_dems<-elections3$'100'
elections3$v31_reps<-elections3$'200'

grouped <- group_by(elections2, v03, v05,v07, v09)
elections2b<-dplyr::summarise(grouped, v30=sum(v30_dems, na.rm=T), v31=sum(v31_reps, na.rm=T))
elections2<-merge(elections2, elections2b, by=c("v03", "v05","v07","v09"))
elections2$v56<-as.numeric(as.vector(elections2$v56))

grouped <- group_by(elections3, v03, v05,v07, v09)
elections3b<-dplyr::summarise(grouped, v30=sum(v30_dems, na.rm=T), v31=sum(v31_reps, na.rm=T))
elections3<-merge(elections3, elections3b, by=c("v03", "v05","v07","v09"))
elections3$v56<-as.numeric(as.vector(elections3$v09))
elections3$v08<-as.character(elections3$v08)
elections3$v16<-as.character(elections3$v16)
elections3$v43<-as.character(elections3$v43)
elections3$v44<-as.character(elections3$v44)
elections3$v45<-as.character(elections3$v45)
elections3$v46<-as.character(elections3$v46)
elections3$v47<-as.character(elections3$v47)
elections3$v48<-as.character(elections3$v48)
elections3$v49<-as.character(elections3$v49)
elections3$v50<-as.character(elections3$v50)
elections3$v51<-as.character(elections3$v51)
elections3$v52<-as.character(elections3$v52)
elections3$v53<-as.character(elections3$v53)
elections3$v54<-as.character(elections3$v54)
elections<-full_join(elections1, elections2)
elections<-full_join(elections, elections3)


elections<-elections[,c(3,6,8,16,17,26,29,34:36)]
elections$demshare<-elections$v30/(elections$v30+elections$v31)
elections$margin<-(elections$demshare-.5)*100
elections$margin_abs<-abs(elections$demshare-.5)*100
colnames(elections)<-c("abb", "election_year","chamber", "district","multimember","party","winner", "totalvotes", "demvotes", "repvotes", "demshare", "margin", "margin_abs")
elections2<-elections[order(elections$margin_abs),]
elections2<-subset(elections2, winner==1)
elections2<-subset(elections2, chamber==9)
elections2$margin_votes<-(elections2$demvotes-elections2$repvotes)/2

elections2$multimember[elections2$abb=="VA" & elections2$year==2014]<-0
elections2$multimember[elections2$abb=="NJ" & elections2$year==2014]<-1
elections2$multimember[elections2$abb=="VA" & elections2$year==2013]<-0
elections2$multimember[elections2$abb=="NJ" & elections2$year==2013]<-1

elections2$multimember2[elections2$multimember==3 | elections2$multimember==6 | elections2$multimember==9 | elections2$multimember==10]<-1
elections2$multimember2[elections2$multimember==1 | elections2$multimember==2 | elections2$multimember==4 | elections2$multimember==5 | elections2$multimember==7 | elections2$multimember==8]<-0

elections2$year<-elections2$election_year

today<-format(Sys.time(), "%Y%m%d")



md_rdd<-matrix(data=NA, nrow=1, ncol=9)
colnames(md_rdd)<-c("abb", "election_year", "margin_seats_abs","margin_seats","margin_votes", "manhattan_distance", "euclidean_distance", "uniform_swing", "multimember")
md_rdd <- as.data.frame(md_rdd)

which(stateleg_elections$abb=="DE" & stateleg_elections$election_year==1980)
for (i in 1:dim(stateleg_elections)[1]){
	if (!is.na(stateleg_elections$margin_abs[i])){
		temp  <- subset(elections2, abb==stateleg_elections$abb[i] & election_year==stateleg_elections$election_year[i])
		## If democrats control the state house
		if ( stateleg_elections$margin[i]> 0){
			temp<-subset(temp, margin>0)
			total_votes<-sum(temp$totalvotes, na.rm=T)
			temp2 <- temp[1:round2(stateleg_elections$margin_abs[i],0),]
			swing<-max(temp2$margin)
			temp2$margin_squared <- temp2$margin^2
			temp3 <- sum(temp2$margin, na.rm=T)
			temp4 <- sum(temp2$margin_sq, na.rm=T)
			margin_votes<- sum(temp2$margin_votes, na.rm=T)
			margin_votes<- margin_votes/total_votes
		}
		## If republicans control the state house
		if ( stateleg_elections$margin[i]< 0){
			temp  <- subset(elections2, abb==stateleg_elections$abb[i] & election_year==stateleg_elections$election_year[i])
			temp<-subset(temp, margin<0)
			total_votes<-sum(temp$totalvotes, na.rm=T)
			temp2 <- temp[1:round2(stateleg_elections$margin_abs[i],0),]
			swing<-min(temp2$margin)
			temp2$margin_squared <- temp2$margin^2
			temp3 <- sum(temp2$margin, na.rm=T)
			temp4 <- sum(temp2$margin_squared, na.rm=T)
			margin_votes<- sum(temp2$margin_votes, na.rm=T)
			margin_votes<- margin_votes/total_votes
		}	
		## If the state house is tied
		if ( stateleg_elections$margin[i]== 0){
			temp  <- subset(elections2, abb==stateleg_elections$abb[i] & election_year==stateleg_elections$election_year[i])
			temp<-subset(temp, margin<0)
			temp2 <- temp[1:round2(stateleg_elections$margin_abs[i],0),]
			swing<- NA
			temp2$margin_squared <- NA
			temp3 <- NA
			margin_votes<-NA
		}	

		temp5 <- c(stateleg_elections$abb[i],stateleg_elections$election_year[i],stateleg_elections$margin_abs[i], stateleg_elections$margin[i], margin_votes,temp3, sum(temp4, na.rm=T)^.5,swing, max(temp$multimember2))
		if (temp3<0 & !is.na(temp3)){
				temp5[7]<-as.numeric(as.vector(temp5[7]))*-1
			}
		if (is.na(temp3)){
				temp5[7]<-NA
			}
		md_rdd <- rbind(md_rdd, temp5)
	}
	print(i)
	if (i==200 | i==400 | i==600 | i==800){
		md_rdd$margin_seats_abs <- as.numeric(as.vector(md_rdd$margin_seats_abs))
		md_rdd$margin_seats <- as.numeric(as.vector(md_rdd$margin_seats))
		md_rdd$manhattan_distance <- as.numeric(as.vector(md_rdd$manhattan_distance))
		md_rdd$euclidean_distance <- as.numeric(as.vector(md_rdd$euclidean_distance))
		md_rdd$uniform_swing <- as.numeric(as.vector(md_rdd$uniform_swing))
		md_rdd$margin_votes <- as.numeric(as.vector(md_rdd$margin_votes))
		
		print(cor(md_rdd$euclidean_distance, md_rdd$margin_seats, use="complete.obs"))
		print(cor(md_rdd$uniform_swing, md_rdd$margin_seats, use="complete.obs"))
		print(cor(md_rdd$margin_votes, md_rdd$margin_seats, use="complete.obs"))
	}
}
dim(md_rdd)
md_rdd$margin_seats_abs <- as.numeric(as.vector(md_rdd$margin_seats_abs))
md_rdd$margin_seats <- as.numeric(as.vector(md_rdd$margin_seats))
md_rdd$manhattan_distance <- as.numeric(as.vector(md_rdd$manhattan_distance))
md_rdd$euclidean_distance <- as.numeric(as.vector(md_rdd$euclidean_distance))
md_rdd$uniform_swing <- as.numeric(as.vector(md_rdd$uniform_swing))
md_rdd$margin_votes <- as.numeric(as.vector(md_rdd$margin_votes))

#md_rdd <- subset(md_rdd, election_year<2011)
#md_rdd$md_rdd[md_rdd$margin < 0 & !is.na(md_rdd$margin)] <- md_rdd$md_rdd[md_rdd$margin < 0 & !is.na(md_rdd$margin)] * -1

cor(md_rdd$euclidean_distance, md_rdd$margin_seats, use="complete.obs")
cor(md_rdd$uniform_swing, md_rdd$margin_seats, use="complete.obs")
cor(md_rdd$margin_votes, md_rdd$margin_seats, use="complete.obs")
md_rdd$year<-md_rdd$election_year

data2<-data
#data2$election_year<-data2$year-1
data2<-dplyr::select(data2,-election_year)
data2<-merge(data2, md_rdd, by=c("year", "abb"), all.x=T)
data2<-subset(data2, year>1966)
#data2<-subset(data2, year<2013)
data2<-subset(data2, hs_dem_per_2pty!=50)
data2<-subset(data2, abb!="DC")
data2<-subset(data2, abb!="NE")
data2<-subset(data2, !is.na(Policy))
#data2<-subset(data2, !is.na(euclidean_distance))


## drop states with MMD
#data2<-subset(data2, abb!="AZ" &abb!="ID" & abb!="MD" & abb!="NH" & abb!="NJ" & abb!="ND" & abb!="SD" & abb!="VT" )

data2<-subset(data2, multimember==0)

cor(md_rdd$manhattan_distance, md_rdd$margin_seats, use="complete.obs")
cor(md_rdd$euclidean_distance, md_rdd$margin_seats, use="complete.obs")
cor(md_rdd$euclidean_distance, md_rdd$manhattan_distance, use="complete.obs")

data2<-dplyr::select(data2, abb, year, euclidean_distance, margin_seats )

#write.dta(data2,file= "statelegrdd_data_160906.dta")
#data2<-read.dta(file= "statelegrdd_data_160906.dta")
## Reload the main data
load('party_control_data_161121.Rdata')
data<-merge(data,data2, by=c("year", "abb"), all.x=T)
################
#### Setup RDD Variables
################

data <- subset(data, !is.na(abb))

data <- plyr:::arrange(data, abb, year)

data <- mutate(data,
             #   RunningVar = running_euclid,
                 RunningVar = euclidean_distance,
                # RunningVar2 = running_euclid * 2,
                RunningVar2 = euclidean_distance * 2,
               TreatmentVar = as.integer(RunningVar > 0),
                GovLib = Policy,
                GovLibL1 = lag(Policy, 1), ## lead 1
                GovLibL2 =  lag(Policy, 2), ## lead 2
                GovLibP1 = lead(Policy, 1), ## lead 1
                GovLibP2 =  lead(Policy, 2), ## lead 2
                GovLibP12 = (GovLibP1 + GovLibP2)/2, ## ave of leads 1-2
                GovLibD1 = GovLibP1 - GovLib, ## delta 1
                GovLibD2 = GovLibP2 - GovLib,  ## delta 2
                GovLibD12 = (GovLibD1 + GovLibD2)/2 ## ave of deltas 1-2
                   )

data <- mutate(data, GovLibL12 = (Policy_L1 + Policy_L2) / 2)
data <- mutate(data, GovLib = Policy)

res.fe <- resid(lm(GovLib ~ abb + factor(year), data=data))
res.l1 <- resid(lm(GovLib ~ Policy_L1, data=data))
res.l2 <- resid(lm(GovLib ~ Policy_L2, data=data))
res.fe.l1 <- resid(lm(GovLib ~ abb + factor(year) + Policy_L1, data=data))
res.fe.l2 <- resid(lm(GovLib ~ abb + factor(year) + Policy_L2, data=data))

data$GovLib.fe0 <- NA
data$GovLib.fe0[as.integer(names(res.fe))] <- res.fe
data$GovLib.fe1 <- NA
data$GovLib.fe1[as.integer(names(res.fe.l1))] <- res.fe.l1
data$GovLib.fe2 <- NA
data$GovLib.fe2[as.integer(names(res.fe.l2))] <- res.fe.l2
data$GovLib.1 <- NA
data$GovLib.1[as.integer(names(res.l1))] <- res.l1
data$GovLib.2 <- NA
data$GovLib.2[as.integer(names(res.l2))] <- res.l2

data <- group_by(data, abb) %>%
    mutate(GovLib.1L1 = lag(GovLib.1, 1), ## lag of first-difference
           GovLib.1L2 = lag(GovLib.1, 2),
           GovLib.2L1 = lag(GovLib.2, 1), 
           GovLib.2L2 = lag(GovLib.2, 2),
           GovLib.fe0L1 = lag(GovLib.fe0, 1),
           GovLib.fe0L2 = lag(GovLib.fe0, 2),
           GovLib.fe1L1 = lag(GovLib.fe1, 1),
           GovLib.fe1L2 = lag(GovLib.fe1, 2),
           GovLib.fe2L1 = lag(GovLib.fe2, 1),
           GovLib.fe2L2 = lag(GovLib.fe2, 2),
           GovLib.1P1 = lead(GovLib.1, 1), ## lead of first-difference
           GovLib.1P2 = lead(GovLib.1, 2),
           GovLib.2P1 = lead(GovLib.2, 1), 
           GovLib.2P2 = lead(GovLib.2, 2),
           GovLib.fe0P1 = lead(GovLib.fe0, 1),
           GovLib.fe0P2 = lead(GovLib.fe0, 2),
           GovLib.fe1P1 = lead(GovLib.fe1, 1),
           GovLib.fe1P2 = lead(GovLib.fe1, 2),
           GovLib.fe2P1 = lead(GovLib.fe2, 1),
          GovLib.fe2P2 = lead(GovLib.fe2, 2),
           GovLibP3 = lead(GovLib, 3),
           GovLibD3 = GovLibP3 - GovLib,  ## delta 3
           GovLibP4 = lead(GovLib, 4),
           GovLibD4 = GovLibP4 - GovLib, 
           DemControlHouse=hs_dem_control,
            DemControlSenate=sen_dem_control,
            HouseControlP4 = (lead(DemControlHouse, 4)) ,
              HouseControlP2 = (lead(DemControlHouse, 2)) ,
          HouseControlP1 = (lead(DemControlHouse, 1)) ,
                HouseControlD2 = (HouseControlP2 -DemControlHouse ) ,
        HouseControlD4 = (HouseControlP4 -DemControlHouse ) ,
             HouseControlD4_1 = (HouseControlP4 -HouseControlP1 ) 

           )


data<-subset(data, !is.na(abb))
data<-subset(data, !is.na(election_year))


################################################################################
#### ESTIMATES #################################################################
################################################################################
####
### MAIN PLACEBOS
#####
with(data, rdrobust(y=gov_party, x=RunningVar2,
                     bwselect="mserd", all=TRUE, kernel='tri'))

with(data, rdrobust(y=DemControlHouse, x=RunningVar2,
                     bwselect="mserd", all=TRUE, kernel='tri'))
##> passed

##> passed

with(data, rdrobust(y=hs_dem_per_2pty, x=RunningVar2,
                    bwselect="mserd", all=TRUE, kernel='tri'))
##> passed

with(data, rdrobust(y=DemControlSenate, x=RunningVar2,
                     bwselect="mserd", all=TRUE, kernel='tri'))
##> passed

with(data, rdrobust(y=sen_dem_per_2pty, x=RunningVar2,
                    bwselect="mserd", all=TRUE, kernel='tri'))
##> passed

with(data, rdrobust(y=Policy, x=RunningVar2,
                    bwselect="mserd", all=TRUE, kernel='tri'))
##> passed

with(data, rdrobust(y=data$GovLib.fe0, x=RunningVar2,
                    bwselect="mserd", all=TRUE, kernel='tri'))
##> passed

with(data, rdrobust(y=GovLib.1L1, x=RunningVar2,
                    bwselect="mserd", all=TRUE, kernel='tri'))
##> passed


#####
## Contemporaneous
######
with(data, rdrobust(y=GovLib, x=RunningVar2,
                     bwselect="mserd", all=TRUE, kernel='tri'))
##> "passed"
with(data, rdrobust(y=GovLib - GovLibL1, x=RunningVar2,
                     bwselect="mserd", all=TRUE, kernel='tri'))
##> passed
with(data, rdrobust(y=GovLib - GovLibL2, x=RunningVar2,
                     bwselect="mserd", all=TRUE, kernel='tri'))
##> "passed"
with(data, rdrobust(y=((GovLib - GovLibL1) + (GovLib - GovLibL2)) / 2,
                     x=RunningVar2, bwselect="mserd", all=TRUE, kernel='tri'))
##> "passed"
with(data, rdrobust(y=GovLib.1, x=RunningVar2,
                     bwselect="mserd", all=TRUE, kernel='tri'))
##> "passed"
with(data, rdrobust(y=GovLib.fe0, x=RunningVar2,
                     bwselect="mserd", all=TRUE, kernel='tri'))
##> passed
with(data, rdrobust(y=GovLib.fe1, x=RunningVar2,
                     bwselect="mserd", all=TRUE, kernel='tri'))
##> "passed"
with(data, rdrobust(y=GovLib.fe2, x=RunningVar2,
                     bwselect="mserd", all=TRUE, kernel='tri'))
##> "passed"

###
## Lag 1
###
with(data, rdrobust(y=GovLibL1, x=RunningVar2,
                     bwselect="mserd", all=TRUE, kernel='tri'))
##> borderline
with(data, rdrobust(y=GovLibL1 - GovLibL2, x=RunningVar2,
                     bwselect="mserd", all=TRUE, kernel='tri'))
##> passed
with(data, rdrobust(y=GovLib.1L1, x=RunningVar2,
                     bwselect="mserd", all=TRUE, kernel='tri'))
##> passed
with(data, rdrobust(y=GovLib.2L1, x=RunningVar2,
                     bwselect="mserd", all=TRUE, kernel='tri'))
##> passed
with(data, rdrobust(y=GovLib.fe0L1, x=RunningVar2,
                     bwselect="mserd", all=TRUE, kernel='tri'))
##> passed
with(data, rdrobust(y=GovLib.fe1L1, x=RunningVar2,
                     bwselect="mserd", all=TRUE, kernel='tri'))
##> passed
with(data, rdrobust(y=GovLib.fe2L1, x=RunningVar2,
                     bwselect="mserd", all=TRUE, kernel='tri'))
##> passed

###
## Lag 2
###
with(data, rdrobust(y=GovLibL2, x=RunningVar2,
                     bwselect="mserd", all=TRUE, kernel='tri'))
##> borderline
with(data, rdrobust(y=GovLib.1L2, x=RunningVar2,
                     bwselect="mserd", all=TRUE, kernel='tri'))
##> passed
with(data, rdrobust(y=GovLib.2L2, x=RunningVar2,
                     bwselect="mserd", all=TRUE, kernel='tri'))
##> passed
with(data, rdrobust(y=GovLib.fe0L2, x=RunningVar2,
                     bwselect="mserd", all=TRUE, kernel='tri'))
##> Borderline
with(data, rdrobust(y=GovLib.fe1L2, x=RunningVar2,
                     bwselect="mserd", all=TRUE, kernel='tri'))
##> margpassedinal
with(data, rdrobust(y=GovLib.fe2L2, x=RunningVar2,
                     bwselect="mserd", all=TRUE, kernel='tri'))
##> passed

##################
### ESTIMATES
##################

#### Main estimates for paper:
## Change t + 1
with(data, rdrobust(y = GovLibD1, x = RunningVar2,
                     bwselect="mserd", all=TRUE, kernel='tri'))
## Change t + 2
with(data, rdrobust(y = GovLibD2, x = RunningVar2,
                     bwselect="mserd", all=TRUE, kernel='tri'))

#### Other estimates:
##> +
## Lead 1
with(data, rdrobust(y=GovLibP1, x=RunningVar2,
                     bwselect="mserd", all=TRUE, kernel='tri'))
##> +
with(data, rdrobust(y=GovLib.1P1, x=RunningVar2,
                     bwselect="mserd", all=TRUE, kernel='tri'))
##> +
with(data, rdrobust(y=GovLib.fe0P1, x=RunningVar2,
                     bwselect="mserd", all=TRUE, kernel='tri'))
##> + **
with(data, rdrobust(y=GovLib.fe1P1, x=RunningVar2,
                     bwselect="mserd", all=TRUE, kernel='tri'))
##> + (similar to FE)
with(data, rdrobust(y=GovLib.fe2P1, x=RunningVar2,
                     bwselect="mserd", all=TRUE, kernel='tri'))

##> + *
## Lead 2
with(data, rdrobust(y=GovLib.fe0P2, x=RunningVar2,
                     bwselect="mserd", all=TRUE, kernel='tri'))
##> + ***
with(data, rdrobust(y=GovLib.2P2, x=RunningVar2,
                     bwselect="mserd", all=TRUE, kernel='tri'))
##> + *
with(data, rdrobust(y=GovLib.fe1P2, x=RunningVar2,
                     bwselect="mserd", all=TRUE, kernel='tri'))
##> + *
with(data, rdrobust(y=GovLib.fe2P2, x=RunningVar2,
                     bwselect="mserd", all=TRUE, kernel='tri'))
##> + ** (similar to FE)

## Average of Lead 1 and 2
with(data, rdrobust(y=(GovLibP1 + GovLibP2) / 2, x=RunningVar2,
                     bwselect="mserd", all=TRUE, kernel='tri'))
##> +
with(data, rdrobust(y=(GovLib.1P1 + GovLib.2P2) / 2, x=RunningVar2,
                     bwselect="mserd", all=TRUE, kernel='tri'))
##> +
with(data, rdrobust(y=(GovLib.fe0P1 + GovLib.fe0P2) / 2, x=RunningVar2,
                     bwselect="mserd", all=TRUE, kernel='tri'))
##> + **
with(data, rdrobust(y=(GovLib.fe1P1 + GovLib.fe1P2) / 2, x=RunningVar2,
                     bwselect="mserd", all=TRUE, kernel='tri'))
##> + **
with(data, rdrobust(y=(GovLib.fe2P1 + GovLib.fe2P2) / 2, x=RunningVar2,
                     bwselect="mserd", all=TRUE, kernel='tri'))
##> +
with(data, rdrobust(y=(GovLib.fe1P1 + GovLib.fe2P2) / 2, x=RunningVar2,
                     bwselect="mserd", all=TRUE, kernel='tri'))
##> + *
## average of change t + 1 and t + 2
with(data, rdrobust(y = GovLibD12, x = RunningVar2,
                     bwselect="mserd", all=TRUE, kernel='tri'))
##> + 

## PRE-1980

with(subset(data, year <= 1980),
            rdrobust(y = GovLibD1, x = RunningVar2,
                     bwselect="mserd", all=TRUE, kernel='tri'))

with(subset(data, year <= 1980),
            rdrobust(y = GovLib.fe1P1, x = RunningVar2,
                     bwselect="mserd", all=TRUE, kernel='tri'))
with(subset(data, year <= 1980),
            rdrobust(y = GovLib.fe2P2, x = RunningVar2,
                     bwselect="mserd", all=TRUE, kernel='tri'))
## POST-1980
with(subset(data, year > 1980),
            rdrobust(y = GovLibD1, x = RunningVar2,
                     bwselect="mserd", all=TRUE, kernel='tri'))

with(subset(data, year > 1980),
            rdrobust(y = GovLib.fe1P1, x = RunningVar2,
                     bwselect="mserd", all=TRUE, kernel='tri'))
with(subset(data, year > 1980),
            rdrobust(y = GovLib.fe2P2, x = RunningVar2,
                     bwselect="mserd", all=TRUE, kernel='tri'))
                     
 ## Graph for paper 

data2<-data
getRD <- function (rdr.out, which.get = c('Conventional', 'Robust')) {
  rdr.out$tabl3.str[which.get, ]
}
myRD <- function (ynm, keep = TRUE, data = data2, xnm = "RunningVar2", c = 0,
                  get = TRUE, ...) {
  rdr <- rdrobust(y = data2[[ynm]], x = data2[[xnm]], subset = keep,
                  c = c, bwselect = "mserd", all = TRUE, kernel = 'tri')
  if (get) {
    getRD(rdr, ...)
  } else {
    rdr
  }
}

n.leads <- 4
(yr.brks <- c(1968, 1989, max(data2$year) - n.leads+2))
data2 <- mutate(data2,
                   Era = cut(year, yr.brks, include.lowest = TRUE, dig.lab = 4))
yr.subs <- llply(levels(data2$Era), function (e) {
  sort(unique(data2$year[data2$Era == e]))
})
names(yr.subs) <- laply(yr.subs, function (yrs) paste0(min(yrs), "-", max(yrs)))
yr.subs[["All Years"]] <- sort(unique(data2$year[!is.na(data2$Era)]))
yr.subs<-list( yr.subs$`All Years`,yr.subs$`1968-1989`, yr.subs$`1990-2012`)
names(yr.subs) <- laply(yr.subs, function (yrs) paste0(min(yrs), "-", max(yrs)))
names(yr.subs)[1]<-"All Years"

(dvs <- paste0("GovLibD", 1:n.leads))
rd.ls <- llply(yr.subs, function (sub) {
  llply(dvs, function (dv) {
    myRD(dv, keep = data2$year %in% sub)
  })
})

rd.cast <- dcast(melt(rd.ls), L1 + L2 ~ Var2 + Var1)
rd.cast <- rbind(rd.cast[9:12,])
rd.cast$L1 <- reorder(rd.cast$L1,1:12)

rd.cast[, -c(1:2)] <- sapply(rd.cast[, -c(1:2)], as.numeric)
rd.cast <- mutate(rd.cast, Lead = factor(L2, labels = 1:n.leads))
summary(rd.cast)

mypdf("Figure5_StateLegRDAll", width=5, height=3)
(ggplot(subset(rd.cast))
 + geom_pointrange(aes(x = Lead, y = Coef_Conventional, ymin = `CI Lower_Robust`,
                       ymax = `CI Upper_Robust`))
 + geom_hline(yintercept = 0, linetype = 'dotted')
 #+ facet_wrap(~L1, ncol = 4)
 + ylab('Effect on Policy Change\n(Robust 95% CI)')
 + xlab("Years after Election")
 + ggtitle("State House")
 + theme_bw()
 )
dev.off()




##> 

(rd.gov <- with(data, rdrobust(y=gov_party, x=RunningVar2,
                                bwselect="mserd", all=TRUE, kernel='tri')))
(rd.hs <- with(data, rdrobust(y=hs_dem_per_2pty, x=RunningVar2,
                               bwselect="mserd", all=TRUE, kernel='tri')))
(rd.sen <- with(data, rdrobust(y=sen_dem_per_2pty, x=RunningVar2,
                                bwselect="mserd", all=TRUE, kernel='tri')))
(rd.L0 <- with(data, rdrobust(y=GovLib, x=RunningVar2,
                               bwselect="mserd", all=TRUE, kernel='tri')))
(rd.fe0 <- with(data, rdrobust(y=GovLib.fe0, x=RunningVar2,
                                  bwselect="mserd", all=TRUE, kernel='tri')))

placebo.mat <- laply(list(rd.gov, rd.hs, rd.sen, rd.L0, rd.fe0), function (x) {
  out <- c(x$coef['Robust', 1], x$se['Robust', 1], x$pv['Robust', 1])
  names(out) <- c('Est', 'SE', 'P>|z|')
  return(out)
})
rownames(placebo.mat) <- c('Dem Governor', 'Dem % in House', 'Dem % in Senate',
                           'Policy Liberalism', 'Residual Policy Liberalism')
xtable(placebo.mat, digits = 3)

## t + 1
(rd.P1 <- with(data, rdrobust(y=GovLibP1, x=RunningVar2, all=TRUE)))
(rd.D1 <- with(data, rdrobust(y=GovLibD1, x=RunningVar2, all=TRUE)))
(rd.P1pre <- with(subset(data, year < 1980),
                  rdrobust(y=GovLibP1, x=RunningVar2, all=TRUE)))
(rd.D1pre <- with(subset(data, year < 1980),
                  rdrobust(y=GovLibD1, x=RunningVar2, all=TRUE)))
(rd.P1post <- with(subset(data, year >= 1980),
                   rdrobust(y=GovLibP1, x=RunningVar2, all=TRUE)))
(rd.D1post <- with(subset(data, year >= 1980),
                   rdrobust(y=GovLibD1, x=RunningVar2, all=TRUE)))
## t + 2
(rd.P2 <- with(data, rdrobust(y=GovLibP2, x=RunningVar2, all=TRUE)))
(rd.D2 <- with(data, rdrobust(y=GovLibD2, x=RunningVar2, all=TRUE)))
(rd.P2pre <- with(subset(data, year < 1980),
                  rdrobust(y=GovLibP2, x=RunningVar2, all=TRUE)))
(rd.D2pre <- with(subset(data, year < 1980),
                  rdrobust(y=GovLibD2, x=RunningVar2, all=TRUE)))
(rd.P2post <- with(subset(data, year >= 1980),
                   rdrobust(y=GovLibP2, x=RunningVar2, all=TRUE)))
(rd.D2post <- with(subset(data, year >= 1980),
                   rdrobust(y=GovLibD2, x=RunningVar2, all=TRUE)))
## ave(t + 1, t + 12)
(rd.P12 <- with(data, rdrobust(y=GovLibP12, x=RunningVar2, all=TRUE)))
(rd.D12 <- with(data, rdrobust(y=GovLibD12, x=RunningVar2, all=TRUE)))
(rd.P12pre <- with(subset(data, year < 1980),
                  rdrobust(y=GovLibP12, x=RunningVar2, all=TRUE)))
(rd.D12pre <- with(subset(data, year < 1980),
                  rdrobust(y=GovLibD12, x=RunningVar2, all=TRUE)))
(rd.P12post <- with(subset(data, year >= 1980),
                   rdrobust(y=GovLibP12, x=RunningVar2, all=TRUE)))
(rd.D12post <- with(subset(data, year >= 1980),
                   rdrobust(y=GovLibD12, x=RunningVar2, all=TRUE)))

est.mat <- laply(list(rd.D12, rd.D12pre, rd.D12post), function (x) {
  out <- c(x$coef['Robust', 1], x$se['Robust', 1], x$pv['Robust', 1])
  names(out) <- c('Est', 'SE', 'P>|z|')
  return(out)
})
rownames(est.mat) <- c('Full Period', '<= 1980', '> 1980')
xtable(est.mat, digits = 3)

summary(group_by(subset(data, year > 1980), year) %>%
            summarise(sd = sd(GovLib, na.rm=TRUE)))

est.mat[3, 1] / 1.1235

################################################################################
#### MAIN EFFECT RD ESTIMATES ##################################################
################################################################################

### RD PLOTS
## Equally spaced
par(mfrow=c(2, 2))
ss <- c(1, 2, 5, 10)
optim.bins <- vector('list', 4)
names(optim.bins) <- paste('times', ss)
for (i in seq_along(ss)) {
  s <- ss[i]
  optim.bins[[i]] <-
    rdplot(y=data2$GovLibD12, x=data2$RunningVar2,
           c = 0, bindplyr::select='es', scale=s, p=4,
           title=paste0('Evenly Spaced Bins, Scale = ', s),
           x.label='Democratic Share t',
           y.label='Change in Policy Liberalism, t to t+2')
  abline(h=0, lty='dotted')
}

## Bins with equal numbers of units (quantiles)
rdplot(y=data2$GovLibD12,
       x=data2$RunningVar2,
       bindplyr::select='qs', scale=1, p=3,
       title='Quantile-Spaced Bins',
       x.label='Democratic Gubernatorial Margin t',
       y.label='Ave. Change in Policy Liberalism, t+1 and t+2',
       y.lim=c(-.11, .11), x.lim=c(-.25, .25))

## Estimates with different "optimal" bandwidths
bw.mserd.tri <- with(data2,
                     rdbwselect(y=GovLibD12, x=RunningVar2, bwselect="mserd",
                                kernel='tri'))
(mserd.tri <- with(data2, rdrobust(y=GovLibD12, x=RunningVar2,
                                      h=bw.mserd.tri$bws[, "h_l"],
                                      b=bw.mserd.tri$bws[, "b_l"],
                                      all=TRUE, kernel='tri')))
## optimal bandwidth for uniform kernel (simple window)
bw.mserd.uni <- with(data2,
                     rdbwselect(y=GovLibD12, x=RunningVar2, bwselect="mserd",
                                kernel='uni'))
(mserd.uni <- with(data2, rdrobust(y=GovLibD12, x=RunningVar2,
                                      h=bw.mserd.uni$bws[, "h_l"],
                                      b=bw.mserd.uni$bws[, "b_l"],
                                      all=TRUE, kernel='uni')))
## OLS with clustered SEs in optimal window (linear)
with(subset(data2, abs(RunningVar2) < bw.mserd.uni$bws[, "h_l"]),
     robcov(ols(GovLibD12 ~ TreatmentVar * RunningVar2, x=TRUE),
            cluster=abb))
## OLS with clustered SEs in optimal window (quadratic)
## According to what Rocio told me, the confidence interval from this should be
## very close to the robust CI for the linear specification.
with(subset(data2, abs(RunningVar2) < bw.mserd.uni$bws[, "h_l"]),
     robcov(ols(GovLibD12 ~ TreatmentVar*pol(RunningVar2, 2), x=TRUE)))
with(subset(data2, abs(RunningVar2) < bw.mserd.uni$bws[, "h_l"]),
     robcov(ols(GovLibD12 ~ TreatmentVar*pol(RunningVar2, 2), x=TRUE),
            cluster=abb))
            

width <- 2
data2 <- mutate(data2, bin=cut(RunningVar2, breaks=seq(-60, 60, width)))

bw.mserd.tri1 <- with(data2,
                       rdbwselect(y=GovLibD1, x=RunningVar2, bwselect="mserd",
                                  kernel='tri'))
(mserd.tri1 <- with(data2, rdrobust(y=GovLibD1, x=RunningVar2,
                                        h=bw.mserd.tri$bws[, "h_l"],
                                        b=bw.mserd.tri$bws[, "b_l"],
                                        all=TRUE, kernel='tri')))

bw.mserd.tri2 <- with(data2,
                       rdbwselect(y=GovLibD2, x=RunningVar2, bwselect="mserd",
                                  kernel='tri'))
(mserd.tri2 <- with(data2, rdrobust(y=GovLibD2, x=RunningVar2,
                                        h=bw.mserd.tri$bws[, "h_l"],
                                        b=bw.mserd.tri$bws[, "b_l"],
                                        all=TRUE, kernel='tri')))

bw.mserd.tri12 <- with(data2,
                       rdbwselect(y=GovLibD12, x=RunningVar2, bwselect="mserd",
                                  kernel='tri'))
(mserd.tri12 <- with(data2, rdrobust(y=GovLibD12, x=RunningVar2,
                                        h=bw.mserd.tri$bws[, "h_l"],
                                        b=bw.mserd.tri$bws[, "b_l"],
                                        all=TRUE, kernel='tri')))

bin1.df <- data.frame(bin.mean=tapply(data2$GovLibD1,
                         data2$bin, mean, na.rm=TRUE),
                    mid=seq(-60 + width/2, 60 - width/2, width),
                     n=tapply(data2$GovLibD1, data2$bin, length))


mypdf("FigureExtra_SH_RDestmserd1-", width=8, height=6)
(ggplot(data2)
 + geom_point(aes(x = RunningVar2, y = GovLibD1), alpha=.1)
 + geom_smooth(data=subset(data2, RunningVar2 > 0),
               aes(x = RunningVar2, y = GovLibD1,
                   weight=tri(RunningVar2, mserd.tri1$h)),
               method = 'lm', formula = y ~ poly(x, 1), size=1.5)
 + geom_smooth(data=subset(data2, RunningVar2 < 0),
               aes(x = RunningVar2, y = GovLibD1,
                   weight=tri(RunningVar2, mserd.tri1$h)),
               method = 'lm', formula = y ~ poly(x, 1), size=1.5)
 + geom_point(data=bin1.df, aes(x=mid, y=bin.mean, size=n), pch=1)
 + geom_vline(xintercept=0)
 + geom_hline(yintercept=0, lty='dotted')
 + coord_cartesian(ylim = c(-.15, .15))
 + scale_x_continuous(breaks=seq(-60, 60, 10), limits=c(-mserd.tri12$h, mserd.tri12$h))
 + labs(size = '# Obs. in\n2-unit Bin')
 + theme_bw()
 + labs(x = "Democrats' Electoral Distance to House Majority Status",
        y = "Change in Policy Liberalism")
 + ggtitle("First Year after Legislative Election")
 + annotate('text', x = 3, y = -.015, 0.05, hjust=0, parse=TRUE,
            label = paste('hat(tau)==', round(mserd.tri1$coef['Conventional', ], 3),
                '~(list(', round(mserd.tri1$ci['Robust', 1], 3),
                ',', round(mserd.tri1$ci['Robust', 2], 3), '))'))
 )
dev.off()




### WINDOWS
yr.window <- 20
window.mids <- seq.int(1968 + yr.window / 2,
                       max(data2$election_year) - yr.window / 2 - 1)

window.ls <- llply(window.mids, function (mid) {
  seq.int(mid - yr.window / 2, mid + yr.window / 2)
})

rd.win.ls <- llply(window.ls, function (win) {
  myRD("GovLibD1", keep = data2$election_year %in% win)
})
names(rd.win.ls) <- window.mids

rd.win.cast <- melt(rd.win.ls) %>% dcast(L1 ~ Var2 + Var1)
rd.win.cast <- sapply(rd.win.cast, as.numeric) %>% as.data.frame()

mypdf("Figure6_StateLegRDWindows-", width=12, height=4)
(ggplot(rd.win.cast)
 + geom_pointrange(aes(x = as.numeric(as.character(L1)),
                       y = Coef_Conventional, ymin = `CI Lower_Robust`,
                       ymax = `CI Upper_Robust`))
 + geom_hline(yintercept = 0, linetype = 'dotted')
 + ylab('Effect on Policy Change\n(Robust 95% CI)')
 + xlab("Year Range")
 + scale_x_continuous(breaks=seq(1940, 2010, 10),
                      labels=paste0(seq(1940, 2010, 10) - 10, "-",
                                    seq(1940, 2010, 10) + 10)) + theme_bw()
 )
dev.off() 


## reverse incumbency advantage

(rd.gov <- with(data, rdrobust(y=HouseControlD2, x=RunningVar2,
                                bwselect="mserd", all=TRUE, kernel='tri')))

